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ABSTRACT 


This thesis presents a field validation of a stochastic, nonlinear wave shoaling 
model based on a third-order closure Boussinesq equations (Herbers and Burton, 1997). 
The model predicts the evolution of directionally spread waves propagating over an 
alongshore uniform beach. The model consists of a coupled set of evolution equations 
for the wave spectrum and bispectrum that incorporates linear shoaling and refraction 
effects and nonlinear energy exchanges in near-resonant triad interactions. Dissipation 
due to breaking is approximated using an empirical quasi-linear damping function and a 
relaxation to Gaussian statistics. The model was verified with field data from five 
alongshore instrument arrays deployed near Duck, North Carolina from August to 
December 1997 as part of the SandyDuck experiment. The predicted shoaling evolution 
of the frequency-directional wave spectra shows the expected development of harmonic 
peaks through triad interactions. The predicted harmonic spectral levels and direction are 
in good agreement with the observed spectra, but the predicted directional spread is 
biased low inside the surf zone. The significant wave height predictions are generally in 
good agreement with observations. The model tends to overshoal the waves outside the 
surf zone and slightly overdissipate wave energy inside the surf zone. Infragravity wave 
growth, sea surface skewness and asymmetry are predicted fairly accurately by the 


model. 
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I. INTRODUCTION 


The focus of this thesis is the verification of a new stochastic Boussinesq wave 
model for the transformation of random surface waves across a natural beach (Herbers 
and Burton, 1997, Herbers et al; manuscript in preparation). Accurate predictions of 
wave shoaling transformation are important for both military (e.g., the planning of 
successful U.S. naval amphibious operations) and civilian (e.g. predicting beach erosion 
during storms) applications. Both linear and nonlinear processes affect the characteristics 
of ocean waves as they shoal onto beaches. Well-established linear theory predicts the 
increase of wave amplitudes, the decrease of wavelengths, and refraction of propagation 
directions toward normal incidence. Changes in directional characteristics of waves are 
particularly important for wave-driven longshore currents and sediment transport in the 
surf zone. In shallow water, nonlinear effects become pronounced, transforming initially 
symmetric, sinusoidal wave profiles to asymmetric, pitched forward profiles of near- 
breaking waves. Near-resonant nonlinear interactions between triads of wave 
components transfer energy from the incident wave components to higher-frequency 
(harmonic) and lower-frequency (infragravity) components. These interactions not only 
broaden the wave spectrum, but also phase-couple the spectral components causing the 
characteristic steepenine and pitched-forward profiles of near-breaking waves (e.g., 
Freilich and Guza, 1984;Elgar and Guza, 1985), and skewed near-bed velocity profiles 
that are believed to play an important role in cross-shore sediment transport. This 
nonlinear evolution process is described well by depth-integrated Boussinesq equations 
for weakly nonlinear, weakly dispersive waves in varying depth (Peregrine, 1967). 


Time-domain Boussinesq models can be applied in principle to natural beaches with 
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arbitrary incident wave conditions, but obtaining complete solutions is computationally 
expensive for large domains and requires detailed boundary conditions that are often 
unavailable. So far, practical applications have been mostly restricted to 1-D 
implementations for uni-directional waves. 

Recently, stochastic models for the transformation of random waves on a beach 
have been developed that solve evolution equations for Statistically averaged spectral 
wave properties (e.g., Abreu et al., 1992; Eldeberky and Battjes, 1996; Agnon and 
Sheremet, 1997; Herbers and Burton, 1997). These models are numerically efficient and 
can be initialized with wave spectra at an offshore boundary, which are easily obtained 
from routine wave measurements or regional wave model predictions. However, 
stochastic models require a statistical closure that may result in significant errors for long 
propagation distances and areas of strong nonlinearity. 

Wheras several theories exist for nonlinear wave interactions, the wave breaking 
process in the surf zone is not yet fully understood and must be described by heuristic 
models. The majority of random wave breaking models are based on the analogy of 
individual wave crests with turbulent bores (Battjes and Janssen, 197 8; Thornton and 
Guza, 1982). These bore models provide robust estimates of bulk dissipation rates within 
the surf zone. However, the spectral characteristics of the energy losses are unknown, 
and somewhat arbitrary quasi-linear spectral forms of the dissipation function are used in 
Boussinesq models (Mase and Kirby, 1992; Eldeberky and Battjes, 1996). 

Despite many studies on the shoaling evolution of wave frequency spectra, the 
evolution of directional wave spectra has not been thoroughly examined. The evolution 


of directional wave spectra is of great importance to the generation of longshore currents 














and sediment transport, and the excitation of nearshore infragravity (0.005-0.05 Hz) 
motions. Laboratory and field experiments of directional wave spectra outside the surf 
zone show the expected refraction of incident waves and energy transfers to harmonic 
components with propagation directions that are in agreement with the theoretical triad 
interaction rules (Freilich et al., 1990; Elgar et al., 1993). Field observations of wave 
transformation across the surf zone show that wave breaking does not significantly alter 
mean propagation directions, but it significantly increases the directional spreading 
(Herbers et al., 1999). The mechanisms for this directional broadening are not currently 
understood. 

This thesis presents a field validation test of a stochastic Boussinesq model 
(Herbers and Biron: 1997) for directionally spread waves propagating over an 
alongshore uniform beach. Norheim et al. (1998) tested a one-dimensional 
implementation of the model for nonbreaking, unidirectional waves with observations of 
wave frequency spectra evolution outside the surf zone on a natural beach. The predicted 
spectra were found to be in good agreement with both measured spectra and predictions 
of the deterministic Boussinesq model of Freilich and Guza (1984). A full two- 
dimensional implementation of the stochastic model for directionally spread waves was 
recently completed (Herbers et al., manuscript in preparation) including a 
parameterization of surf zone dissipation (Whitford, 1988). 

The field data used for model validation were collected during the SandyDuck 
experiment conducted near Duck, North Carolina, during the fall of 1997 (Feddersen et 
al., 2000; Elgar et al., 2001). An extensive array consisting of 69 pressure gages, 33 


electromagnetic current meters, and 33 downward-looking sonar altimeters to monitor the 














sea bed level, was deployed in a 200 m (alongshore) by 420 m (cross-shore) area 
extending from about 1 to 5 m water depth. High quality data were collected during a 
four-month period spanning a wide range of conditions. Four alongshore pressure sensor 
arrays were selected from this experiment to estimate the evolution of the frequency- 
directional wave spectrum E(f,6) across the beach. Additional estimates of E(f,@) were 
extracted from the permanent pressure gage array maintained in 8m depth by the Army 
Corps of Engineers. 

The model was initialized with £(/,6) estimates in both 8m and 5m depth to 
examine the sensitivity of the model to initial conditions. Model predictions of E(f,@) at 
the shallower arrays were compared with the observed spectra. The shoaling evolution of 
several important wave parameters will be examined to quantify the accuracy of the 
model. These parameters include the significant wave height, sea surface skewness and 
asymmetry, infragravity variance, mean propagation direction, and directional spread. 
Several case studies will be examined in detail to better understand model tendencies. 
The overall model performance will be examined using statistical analysis of predictions 
derived from the entire experiment. 

A description of the field experiment and procedures for data reduction and 
analysis are given in chapter II. The model equations and implementation are briefly 
reviewed in chapter III. Several case studies with different wave conditions are examined 
in chapter IV. Results for the entire data set are presented in chapter V, followed by 


conclusions in chapter VI. 














Il. FIELD EXPERIMENT AND DATA ANALYSIS 


The non-linear evolution of the frequency-directional wave spectrum E(/,@) 
across a natural beach is examined using extensive measurements collected at the U.S. 
Army Corps of Engineer’s Field Research Facility (FRF), Duck, NC, during the 
SandyDuck experiment in the summer/fall of 1997 (Elgar et al., 2001; Feddersen et al., 
2000). This site is located on a straight barrier island with an 80 km wide, shallow (20- 
50 m depth) shelf that is fully exposed to the Atlantic Ocean. This study will focus on 
the densely instrumented nearshore region (210-900 m offshore) where the shallow water 
effects on the evolution of the E(/,6) (nonlinear triad interactions, surf zone breaking) 
are most pronounced. 

The bathymetry (figure 1) 1s characterized by a gently sloping inner shelf 
(~1:250) and a slightly steeper beach (~1:100). Bathymetric surveys of the nearshore 
region indicate that the beach topography changed slightly over the course of the 
experiment. The sand bar located approximately 300m offshore remained fairly 
stationary with a crest approximately 5 m below mean sea level. The beach morphology 
was more dynamic on the beach face shoreward of the nearshore arrays, but this region is 
not examined in this study. 

A plan view of the beach with instrument locations is shown in figure 2. 
Directional wave data were available from a permanent array of pressure sensors located 
about 900 m from shore in 8 m depth, which is maintained by the FRF. This 15-element 
array has an aperture of 250 m alongshore and 120 m cross-shore, and is sampled at a 2 
Hz rate that resolves waves with frequencies up to about 0.3 Hz. A subset of 9 sensors 


aligned in the alongshore direction was used here to estimate E(/,@) (figure 2). 
; | 

















The extensive 2-dimensional nearshore array consisted of 69 pressure gauges (P) 
and 33 current meters (U,V). Four linear alongshore arrays of pressure sensors each with 
an aperture of 210 m were used to obtain estimates of E(/,6). These arrays (labeled B, 


C, E, and F in figure 2) were positioned at offshore distances of 210, 260, 375, 500 m in 
nominal depths of 2.9, 3.7, 3.9, and 5.1 m, respectively (figure 2). The nearshore arrays 
also included downward-looking sonar altimeters (S) that provided continuous sea-bed 
level measurements. All nearshore instruments were connected to a central data 
collection center onshore. Data were collected nearly continuously with a sample rate of 
2 Hz between 2 August and 3 December, 1997. 

Alongshore homogeneity over the extent of the nearshore arrays was verified by 
intercomparing the measured spectra of each sensor in the array, as well as the coherence 
and phase spectra of redundant array lags. Elgar et al. (2001) noted that the southernmost 
cross-shore sensor transect of the array consistently recorded spectral density values 
lower than other array sensors for waves arriving from large, oblique southerly directions. 
These discrepancies are likely caused by the close proximity to the FRF research pier 
pilings creating partial wave blocking and a large scour hole (figure 1). Data from these 
sensors were excluded from the analysis. 

In this study, frequency-directional wave spectra were estimated from five 
alongshore arrays of pressure sensors (figures 2,3). Cross-shore arrays were not used 
because cross-shore depth variations cause Petites: cross-shore variations in spectral | 
levels, particularly in surf zone where strong dissipation causes large gradients in wave 
height. The 180 degree directional ambiguity of the linear arrays was resolved by 


neglecting weak reflection of waves from shore (Elgar et al., 1994). The data from the 




















arrays were archived in synchronized 3-hour data records. Only the first 102 minutes of 
each record were used here to estimate E(f,@) 1n order to reduce nonstationary effects of 
tidal fluctuations in water depths. 

Array cross-spectra were calculated based on Fourier transforms of overlapping 
512 sec segments. Merging 3 frequency bands resulted in estimates with a frequency 
resolution of 0.0059 Hz and 72 degrees of freedom. At each frequency, /, the 
directional distribution of wave energy, S(@;), was estimated from the array cross- 
spectra using the variational technique described in Herbers and Guza (1990). The 
direction, @ , is defined relative to the local shoreline Sapaeien with @ =0 
corresponding to normal incidence to the beach (waves arriving from 70 true N), and @ is 
positive (negative) for waves approaching the beach from northerly (southerly) 
directions. The surface elevation frequency spectrum, E(/), was obtained by averaging 
the measured auto-spectra and applying a linear theory depth correction. The S(6; f) 
estimates were combined with E(/) estimates to form the wave frequency-directional | 
spectrum, E( f,0)=E(f)S(6;f). Additionally, bispectra were calculated from the time 
series data ‘s obtain estimates of the bulk wave skewness and asymmetry (e.g. Elgar and 
Guza, 1985). 

Bottom profiles used in numerical model computations were estimated based on 
both the Coastal Research Amphibious Buggy (CRAB) beach surveys conducted by the 
FRF staff and the sonar altimeter data. CRAB high-spatial resolution surveys were 
conducted intermittently in time with gaps ranging from several days to a week. All 
depths were referenced to the NGVD datum. The cross-shore profile in the array vicinity 


measured by the CRAB was linearly interpolated over time to correspond to the spectral 
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wave data estimated at three-hour intervals. Since the CRAB surveys cover only the 
nearshore region out to about 5 m depth, the profiles were extended seaward to the 8 m 
array using a linear interpolation, assuming that the seabed level at the center of the 8 m 
array (-8.08 m) did not change during the course of the experiment. Tide corrections 
were applied to the beach profiles using a 102 minute sea level average obtained from the 
FRF tide gauge. The FRF tide gauge did not operate during the periods 10/21/97 19:00 
EST — 10/22/97 16:00 EST and 11/23/97 19:00 EST — 11/25/97 16:00 EST. The 8m 
array pressure sensor data was used to fill in the missing tide corrections. 

The sonar altimeters recorded sea-bed level data nearly continuously in time 
(same sampling scheme as the pressure sensor data). Only 18 three-hour records were 
missing in the altimeter data during the entire experiment. However, the cross-shore 
resolution was poor with bed-level measurements at only 8 cross-shore positions (see 
figure 4). At each cross-shore position, only the working sonar altimeter located nearest 
the central axis of the arrays (alongshore coordinate 830 m) was used in the analysis. 
Sonar altimeter-based depth profiles were interpolated both spatially and temporally (to 
cover gaps in biestiata collection) and tide corrected. A linear interpolation was used in 
space and time, with the exception of the spatial interpolation across the bar region where 
a cubic spline interpolation was employed to better represent the sand bar crest (figure Sa, 
Sb). 

The root-mean-square (RMS) error between beach profiles estimated from the 
CRAB and altimeter data are shown in figure 6b. The RMS difference between the two 
nearshore profile data sets is generally small (0.04 — 0.08 m), with the exception of a few 


periods (notably 15-23 October and 6-9 and 13-15 November) with more energetic waves 











when the RMS differences were as large as 0.16 m (figures 6a and 6b). In these storm 
events when the CRAB could not safely take measurements, systematic differences 
between the CRAB and altimeter depth estimates suggest significant erosion close to 
shore. In the present study, the CRAB surveys were used in low-moderate wave 
conditions (Hs < 2 m) when sea bed levels did not vary significantly between surveys. 
Profiles based on the sonar altimeter data were used during two wave events (10/18/97 
0100 — 10/20/97 2200; 11/01/97 0100 — 1000) when significant wave heights exceeded 2 


m and the bottom morphology was more dynamic. 
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It. MODEL DESCRIPTION AND IMPLEMENTATION 


The nonlinear, spectral wave shoaling model used in this study is a stochastic 
Boussinesg model for directionally spread waves shoaling on a gently sloping beach with 
straight and parallel depth contours (Herbers and Burton,.1997; Herbers et al., manuscript 
in preparation). The model is based on a third-order closure scheme and consists of a 


coupled set of first-order evolution equations for the wave spectrum E(@,/) and 
bispectrum B(@,,/,,@,,/,) 


dE (ol) _{ 1 dh 
re Oh de D,(@, Ele l) 


(3.1) 


* am 7 [ao far IM {B(o',.@-o@',1-1} 


aB(o'",l',o-o',1-1') 


ax 


{| Re o@'(@ - @'\o _ gh)" (al'- a'ly 


= 328 eR: (a',l',@-o',1-1')- = B(o',l',a@-o',1-1’) 
2¢7" 20'(@- @')@ 


4h dx 
i sseqelo' Elo 0'1-1E(0,))+(@-0B(0'N)E(0,1)- 0F(0'N)E(0~ 01-19) 
(3.2) 

where @ is frequency, / is alongshore wavenumber, / is water depth, g is gravity, and 
IM { } is the imaginary part. The two-dimensional spectrum E(@,/) describes the surface 
elevation spectral density of component (a,/) , and the four-dimensional (complex) 
bispectrum B(,,/,,@,,/,) defines the phase relationship of a triad consisting of 
components (,,/,), (@,,/,), and (-@, -@,,-/,-1,). A heuristic parameterization of wave 
breaking effects was included in equations 3.1 and 3.2 to extend the model into the surf 


zone. The damping terms, D, and D, are expressed in terms of a frequency-dependent 


dissipation function G(@). 
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D,(,1) = G(@) + G(-a) (3.3) 

D,(@',l',@-@',1-1') = G(@') + G(@-0') + G(-a) (3.4) 

Following Mase and Kirby (1992) and Chen et al. (1997) a quadratic frequency 
dependence was included in the model 


G(o,)l) ao’ (3.5) 


The bulk dissipation rate, ¢ = 2pg**h'” [ G(w)dw , was estimated using a parameterization 
p p 


—-® 


given by Whitford (1988) 


We petit, 1am {Ze eG (3.6) 


where f is the mean frequency, H.. is the root mean Square wave height and 7, 5 are 
adjustable coefficients. Additionally, a relaxation term with an adjustable coefficient R is 
introduced in the bispectrum equation to provide a rerandomization of phases in the wave 
breaking process. 

The adjustable breaking parameters were calibrated using extensive observations 
of the evolution of wave frequency spectra across a barred beach (neglecting directional 
effects), and fixed at values: y = (0.30), 6=(0.25), and R= (2.5) (Herbers et.al; manuscript 
in preparation). 

The model is initialized at the offshore boundary with an observed frequency- 
directional spectrum £(/,6) that is transformed to (a, !) space. The initial bispectrum is 
estimated with second order, finite depth theory (Hasselmann et al.,1963). To examine 
the sensitivity to offshore boundary conditions, the model was initialized at both the 8 m | 


and F arrays. 
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Equations 3.1 and 3.2 were integrated using a fourth order Runge Kutta scheme 
with a fixed step size of 0.25 m from the offshore boundary to the B array. The model 
was implemented with 54 frequencies (spanning the full range 0.00585 — 0.3164 Hz of 
the observed spectra), and 85 alongshore wavenumbers. To prevent instability of the 
model, the alongshore wavenumber was restricted (maximum / values of +0.19 rad/m). 
Asa ait of this truncation, high frequency waves entering the model domain at large 
oblique angles are excluded. 

The frequency-directional spectra predicted at all shoreward arrays are compared 
to observed spectra in chapter 4. To quantify the model performance, various important 
bulk parameters: the significant wave height (Hs ), skewness and asymmetry (e.g. Elgar 
and Guza, 1995), infragravity variance (in the frequency range 0.00585 — 0.0468 Hz), and 
the mean direction (@ ) and directional spread («, ) (e.g. Herbers et al., 1999), were 
extracted from the predicted spectra and bispectra, and are compared to observed values 


in chapters 4 and 5. 
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IV. CASE STUDIES 


Four case studies spanning a wide range of wave conditions were selected to 
illustrate in detail the model’s characteristics. In each case, results are presented for 
model runs initialized in both 8m and 5m (F-Array) depths. Case I (August 7, 1997 
01:00 EST) is characterized by a broad spectrum of small, non-breaking waves. Case II 
(August 10, 1997 04:00 EST) features the evolution of a low energy, narrow swell 
spectrum again with little or no breaking. In case III (September 12, 1997 07:00 EST) 
more energetic swell was observed with waves breaking across the inner part of the 
nearshore array transect. Finally in case [TV (October 19, 1997 16:00 EST), collected 
during a severe Noreaster storm that produced the largest waves of the entire data set, the 
entire instrumented transect was within the surf zone. For each case study, the frequency 
spectrum E(f) , and frequency-directional spectrum E(/,6) predicted at each of the 
arrays are compared to the observed spectra (figures 7, 10, 13, 16; only model runs 
initialized at the F array are shown). The predicted and observed cross-shore evolution of 
the significant wave height Hs , sea surface asymmetry and skewness, infragravity 


variance, mean direction 6, and directional spread o, are compared in figures 8, 9, 11, 
12, 14, 15, 17, and 18. The significant wave height (Hs =4,] [fz f, eafae} represents the 


average height of the highest 1/3 waves. Skewness and asymmetry characterize the 
nonlinearity of the waves. The asymmetry describes the pitched forward (positive) or 
backward (negative) profile of the wave, while the skewness describes the peakedness of 
the wave (see Elgar and Guza, 1985 for definitions). The infragravity variance is the 


variance within the infragravity frequency range (0.0059 — 0.0468 Hz) that contains low 
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frequency waves excited by nonlinear interactions of the incident waves in shallow water. 
The mean direction and directional spread are energy weighted bulk average values over 
the spectrum (see Herbers et al., 1999 for definitions). 


CASE I: AUGUST 7, 1997 01:00 EST 


The small incident waves (Hs = 0.52m) are characterized by a bimodal spectrum that 
consists of a southeast swell (0.13 Hz) and a northeasterly swell (0.2 Hz). The predicted 
frequency-directional spectra show virtually no evolution across the transect, in good 
agreement with the observed spectra (figure 7). 

The significant wave heights are accurately predicted in the F array-initialized 
model run, whereas the 8m-initialized run overpredicts Hs by about 10% at the nearshore 
arrays, probably because the high frequency incident waves violate the shallow water 
approximation of the Boussinesq model (figure 8). Predicted sea surface asymmetry and 
skewness values are small (<0.1) at all arrays, consistent with the observed values (figure 
8). 

The observed amplification of in astatiey variance near the shore (about a factor 
of three between the 8m and B arrays) is predicted well by both model runs (figure 9). 
The 8m-initialized model run is slightly more accurate than the F array model run, but 
both model runs underpredict the sharp increase at the B array, possibly because 
reflection from shore is neglected in the model. The observed and predicted mean 
directions are close to normal incidence due to averaging over 2 swell systems arriving 
from opposite quadrants (figure 9). Both model runs reproduce the slight decrease in 


directional spread towards the shore caused by refraction (figure 9). 
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CASE II: AUGUST 10, 1997 04:00 EST 
This day was characterized by moderate size non-breaking waves (Hs = 0.87 m 


offshore increasing to 0.97 m inshore) arriving at near normal incidence (6 = —10 deg). 
The frequency spectrum is narrow with a peak frequency fp = 0.09 Hz (figure 10). 
Strong nonlinear energy transfers are evident in the dramatic growth of harmonic peaks at 
frequencies 2fp = 0.18 Hz and 3fp = 0.27 Hz between the F and C arrays. The predicted 
harmonic spectral levels and directions (aligned with the dominant swell) are in good 
agreement with the observations (figure 10). Interestingly, both the observed and 
predicted spectra show a decay of the second harmonic peak between the C and B arrays. 

Both model runs indicate a slight overshoaling of significant wave heights (figure 
11). The errors are larger for the 8m-initialized model run (about 10% at array B) than 
for the F array-initialized model run (7%), suggesting the errors are caused by 
a of the shallow water approximation used in the Boussinesq equations. 
Asymmetry 1s predicted well by both model runs (figure 11), showing a positive 
maximum of 0.3 (pitched forward wave crests) seaward of the bar and a negative 
maximum of -0.2 (pitched backward crests) inshore of the bar. The observed and 
predicted negative asymmetry indicates that the observed decay of the second harmonic 
peak inshore of the bar is the result of a reversal of nonlinear energy transfers (Norheim 
et al., 1998). Both model runs overpredicted skewness by about 25%, but reproduce 
qualitatively the increase to a maximum value of about 1 over the bar followed by a 
decrease towards the shore (figure 11). 

Again, the infragravity variance increases in the observations as well as the 
models (figure 12). As in case I, the 8m model run predicts more accurately the 


amplification of infragravity waves near the shore, but both model runs underpredict the 
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amplification of infragravity waves at the B array. The predicted shifting of 6 from —10 
deg at the 8m array to —5 deg at the B array is in good agreement with the observed 


directions (figure 12). However, whereas the model predicts a decrease of o, from 15 
deg to 10 deg, the observed o, are nearly uniform across the instrument transect. 


CASE Il; SEPTEMBER 12, 1997 07:00 EST 


The incident waves on this day are characterized by a similar narrow swell 
spectrum (fp = 0.08 Hz, 6= -10 deg), but are more energetic than in case II. The 
significant wave height Hs decreases from a maximum value of 1.4 m offshore of the bar 
to 1.2 m inshore of the bar (figure 14), indicating breaking of the larger waves on the Bar. 
The growth of the harmonic peaks through nonlinear wave-wave interactions is again 
evident in the spectra, and the spectral levels and direction (aligned with the dominant 
swell direction) are well predicted by the model. 

As in the earlier cases, the model tends to overshoal the waves, producing 
maximum wave heights at the E array that are slightly larger than observed. The errors 
are again larger for the model run initialized at the 8m array (8%) than for the run 
initialized at the F-array (5%), suggesting the errors result from inaccuracies in the 
shallow water approximation. The decrease of the significant wave height inshore of the 
E array is accurately predicted by both model runs (figure 14). Predicted asymmetry 
increases to a maximum positive value of 0.4 on the bar (associated with energy transfers 
to high frequency harmonics) followed by a decrease to small negative values inshore of 
the bar (figure 14). The skewness predictions show two maxima (values about 1 at cross- 
shore distances of 250 and 410 m) with slightly smaller values (0.8) on the bar (figure 
14). Predicted skewness and asymmetry variations in both model runs are in good 


agreement with the observations. 
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The amplification of infragravity waves near the shore is predicted accurately by 
the model run initialized in 8m depth, whereas the run initialized at the F array 
underpredicts the infragravity variance by about 30% (figure 15). A possible explanation 
for this large discrepancy is inaccuracy of the initial bispectrum that is based on a local 
horizontal bottom solution of second order wave theory (Herbers and Burton, 1997). 
Both model runs overpredict changes in 6 toward normal incidence by 3-5 degrees and 
overpredict the directional narrowing of the spectrum. The observed co, values at the B 


array are about 30% larger than the predicted values, suggesting that wave breaking 
causes a directional broadening that is not represented in the model. 


CASE IV: OCTOBER 19, 1997 16:00 EST 


In this case with the largest waves of the entire data set, the entire instrument 

_ transect was in the surf zone. The observed and predicted frequency-directional spectra 
again show development of harmonic peaks aligned with the dominant 0.10 Hz swell. 
However, at the inshore C and B arrays predicted nonlinear interactions have filled in the 
valleys between the harmonic peaks, producing a broad, featureless spectrum, in good 
agreement with the observed frequency spectra (figure 16). 

The significant wave height decreases from 3.5 m at the 8m array to 1.5 m at the 
| B array (figure 17). The wave height decay is predicted well by both model runs (figure 
17). The 8m-initialized model run initially predicts too much dissipation 
(underpredicting Hs by 25% at the F array), but yields accurate Hs predictions at the 
inshore arrays (figure 17). The asymmetry and skewness are predicted well by both 
model runs (figure 17). In contrast to nonbreaking case II, the observed and predicted 
asymmetry values are positive over the entire transect, indicating that no reversal in 


energy transfers to lower frequencies occurred in this case. Both model runs predict a 
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maximum asymmetry value at about 0.4 near the E array decreasing to near-zero values 
at the B array, in good agreement with the observations (figure 17). Both model runs 
show nearly uniform (about +0.5) skewness values throughout the transect that are 
slightly (10-20%) smaller than the observed values (figure 17). 

The evolution of the infragravity variance is predicted well by both model runs 
(figure 18). The observed and predicted cross-shore variations are weak compared to the 
previous cases. The infragravity variance decreases near the shore, possibly due to some 
dissipation (figure 18). The predicted turning of 9 toward normal incidence (from +10 
deg in 8m depth to +4 deg at the B array) owing to refraction is in good agreement with 
the observed directions (figure 18). Whereas the predicted o, decreases from 18 degrees 
in 8m depth to 15 degrees at the B array, the observed o, increase slightly to a maximum 
value of 22 degrees at the B array (figure 18). These results suggest that wave breaking 


causes a significant directional broadening of wave spectra in the surf zone. 
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VV. OVERALL MODEL PERFORMANCE 


The overall performance of the model for the entire data set is examined here with | 
scatterplots of predicted versus observed significant wave heights (figure 19), 
infragravity variance (figure 20), skewness and asymmetry(figure 21), and mean 
direction and directional spread (figure 22). 

First, all 987 data runs were examined for bad data. The entire data run was 
discarded if the initialization spectrum from either the 8m or F arrays was missing or bad 
data. A total of 54 data runs were discarded. If data were missing or bad at any one of 
the shoreward E, C, or B arrays, the observed spectrum for that array was discarded. A 
total of 41 spectra were discarded. After bad data were removed, scatterplots were 
created for each array located shoreward of the array used to initialize the model. In each 
scatterplot, a best fit line obtained through linear regression analysis is indicated (dashed 
diagonal line) as well as the one-to-one correspondence line (solid diagonal line). Results 
are shown only for model runs initialized at the F array that yielded the best overall wave 
height predictions (discussed below), with the exception of infragravity variances for 
which both 8m and F array initialized results are shown in figure 20. 

Significant wave height scatterplots (figure 19) show generally good agreement 
for the entire data set, results for both the absolute Hs and Hs normalized by the initial 
Hs are shown. However, the model tends to underpredict Hs (by as much as 10%) when 
the observed Hs is greater than 1.5 m, suggesting a deficiency of the dissipation 
parameterization in high energy conditions. The normalized scatterplots show a tendency 
of a slight overprediction of Hs for nonbreaking waves (normalized Hs >1) and a slight 
underprediction of Hs within the surf zone (normalized Hs <1) (figure 19). 
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Comparisons of observed and predicted (figure 20, both normalized by the intial 
values) infragravity variance show large discrepancies at the arrays located farther 
offshore, but improved correlation as the model progresses shoreward. These results 
suggest the main source of errors is the initialization of the bispectrum that is based on 
second-order wave theory for a horizontal sea bed. The model run initialized at the 8m 
array yields more accurate predictions of the infragravity variance at all shoreward arrays 
than the model run initialized at the F array, consistent with smaller bottom slope effects 
in deeper water. 

The skewness predictions are well correlated with the observations for the entire 
data set, with slightly more scatter at the B array than at the other two arrays (figure 21). 
However, deviations of the best fit line from a one-to-one correspondence show a 
consistent 10-15% overprediction of skewness at all array locations. Skewness (observed 
and predicted) are positive at all arrays indicating that waves always have peaked crest 
profiles in the nearshore region. The asymmetry predictions are also in good agreement 
with observations. The aeyannieey comparisons show slightly more scatter than the 
skewness comparisons, but virtually no bias. Only at the shallower B array, asymmetry 


predictions are biased low for large asymmetry values. Asymmetry is consistently 





positive at the E array, but negative values are observed and predicted at C and B arrays 
(figure 21) that are indicative of reversals of energy transfers to the spectral peak inshore 
of the sand bar. 

Mean direction predictions are well correlated with the observations at all array 
locations (figure 22). The predicted angles are on average slightly closer to normal 


incidence, possibly because high-frequency waves traveling at large oblique angles are 
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excluded in the initial model spectra. Wave refraction is evident in the smaller range of 
observed and predicted 6 at shallow inshore array locations. The directional spread o, 
is consistently underpredicted at all locations (figure 22). This bias, increasing from 
roughly 10% at array E to 25% at array B, is likely caused by the exclusion of high- 
frequency, high angle waves in the model and surf zone scattering effects not represented 


in the model. 
To summarize the accuracy of the model, errors in the model predictions are 


quantified with a mean bias, RMS error, and scatter index (SI), defined as 


N 
bias = +S y, -x, (5.1) 
N ia 


RMS ' error = 





yee! —x;,)’ (5.2) 


_ RMS error 


[L$ 
N i=l 


SI (5.3) 


where y is predicted, x is observed, and N is the number of observations. These error 
statistics were calculated for significant wave height, asymmetry, skewness, infragravity 
variance, mean direction, and directional spread, based on predictions at all arrays 
shoreward of the array used to initialize the model. The results together with a 
correlation coefficient, are listed in tables 1 (model initialized at the 8m array) and 2 


(model initialized at the F array). 
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The F array initialized model runs yielded consistently more accurate significant 
wave heights than the 8m array initialized model runs. This can be deduced from the 
higher correlation coefficient (0.991 for F, 0.986 for 8m) and the lower bias (0.021 m for 
F, 0.047 m for 8m), RMS error (0.055 m for F, 0.081 for 8m), and SI (0.063 for F array, 
0.092 for 8m). Possible explanations for these differences include inaccuracies in the 
shallow water approximation that affect model results in deeper water and a breakdown 
of the dissipation parameterization in large wave conditions. 

Larger errors are observed in sea surface skewness and asymmetry predictions, 
but no significant differences are noted between model runs initialized at the 8m and F 
arrays. 

Infragravity variance is predicted more accurately by the model run initialized in 
8m depth than the model run initialized at the F array as noted in the smaller negative 
bias (-0.00018 m” for 8m, -0.00034 m? for F) , RMS error (0.00056 m2 for 8m, 0.00070 
m2 for F), and SI (0.19 for 8m, 0.21 for F). These differences are likely caused by errors 
in the initialization of the bispectrum. 

Mean direction predictions of model runs initialized at the F atray are more 
accurate than those initialized in 8m depth as exemplified by its higher correlation 
coefficient (0.972 for F, 0.915 for 8m), and the lower RMS error (3.8 deg for F, 5.9 for 
8m) and SI (0.28 for F, 0.40 for 8m). Note that the bias is small owing to cancellation of 
a negative bias for negative (southerly) angles and a positive bias for positive (northerly) 
directions (figure 22). 

There is not much difference between the error statistics of 8m-initialized and F 


array-initialized model predictions of o,. The bias is slightly smaller for the F array- 
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initialized model runs, but in both cases low correlation coefficients (0.792 and 0.807) are 
observed, suggesting an important physical mechanism is missing in the model. Whereas 
the skewness predictions show smaller SI (0.34-0.37) and higher correlation coefficients 
(0.947-0.949) than the asymmetry predictions (0.42-0.45, 0.903-0.905), the asymmetry 
predictions have smaller bias. An important source of errors for both skewness and 
asymmetry predictions is the parameterization of wave breaking. Earlier studies (Chen et 
al., 1997) using a deterministic Boussinesq model show comparable agreement of 
skewness and asymmetry predictions with field and lab data and also indicate a strong 


sensitivity to the choice of dissipation function. 
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VI. SUMMARY 


A stochastic, nonlinear wave shoaling model based on Boussinesq equations and a 
third-order, weakly non-Gaussian closure scheme (Herbers and Burton, 1997) was 
validated using a comprehensive field data set. The model predicts the evolution of 
directionally spread waves propagating over an alongshore uniform beach. The model 
consists of a coupled set of evolution equations for the wave spectrum and bispectrum 
that incorporates linear shoaling and refraction effects and nonlinear energy exchanges in 
near-resonant triad interactions. Dissipation due to breaking is approximated using a 
frequency dependant empirical quasi-linear damping function and a relaxation to 
Gaussian statistics. The bulk dissipation rate is based on a bore model (Thornton and 
Guza, 1983; Whitford, 1988). The model was verified with field data from five - 
_alongshore arrays deployed at the U.S. Army Engineer Field Research Facility, Duck, 
North Carolina from August to December 1997 as part of the SandyDuck experiment. 
Estimates of frequency-directional spectra were obtained from all 5 arrays at 3 hour 
intervals for the entire experiment. Accurate depth profiles were obtained from frequent 
surveys and continuous sea bed level measurements at the array locations. Model runs 
are presented initialized in both 8m or 5m depths to examine the sensitivity of the model 
to initial conditions. 

The shoaling evolution of the frequency-directional spectra was examined in 
detail for four case studies including nonbreaking and breaking waves. Predictions of 
important bulk parameters, significant wave height, sea surface asymmetry and skewness, 
infragravity variance, mean direction, and directional spreading, are compared with 


observed values for the entire experiment. 
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Predicted frequency-directional wave spectra reproduced accurately the observed 
development of harmonic peaks aligned with the dominant wave direction. However, the 
evolution inside the surf zone indicates that the observed spectra are directionally broader 
than the predicted spectra. 

Significant wave height predictions are in good agreement with observations. The 
scatter index is 0.092 for model runs initialized in 8m depth and 0.063 for model runs 
initialized in 5m depth. The comparisons show a systematic model tendency to slightly 
overshoal the waves outside the surf zone (particularly for model runs initialized at the 
deeper 8m array) (tables 1,2). The model also tends to slightly over dissipate the waves 
inside the surf zone. 

Infragravity waves are low frequency waves (nominally 0.005-0.05 Hz) that are 
believed to be important to sand bar formation and beach erosion. The observed dramatic 
infragravity wave growth in shallow water is reasonably well predicted by the model 
except near the shoreline where the model is biased low, possibly owing to neglected 
reflection. Model runs initialized in 8m depth are consistently more accurate than those 
initialized in Sm depth, indicating a sensitivity to the initial conditions that are based on 
second-order wave theory for a horizontal sea bed 

Sea surface skewness and asymmetry are reasonably well predicted by the model, 
but the scatter index values are notably higher than for the wave height predictions 
(skewness 0.37-0.37; asymmetry 0.43-0.45) (tables 1,2). Both observed and predicted 
skewness are consistently positive and the model is biased high. The observed and 
predicted asymmetry are positive (pitched forward wave crests) at the offshore 


measurement locations. Negative asymmetry values (pitched backward wave crests) are 
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often observed and predicted inshore of the sand bar indicating reversal of nonlinear 
energy transfers from harmonics back to the dominant waves. 

Although the stochastic Boussinesq wave model generally performed well, further 
model validation is needed on different beaches. Further possible extensions and 
refinements of the model include allowing for alongshore depth variations and reflection 
of low frequency waves from the shore. As faster computers with more memory become 
available in the future, this efficient spectral model may become a useful tool for real 


time surf forecasts. 
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Figure 1. Bathymetry at the field site surveyed on September 16, 1997. The 
instrumented region is indicated with a red box. The trench located to the south of the 
instrumented region is caused by the presence of a 600m-long pier. 
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Figure 2. Plan View of Arrays. Instrument types include SPUV (collocated sonar, 

pressure sensor, and current meter), PUV (collocated pressure sensor and current meter) 

and P (pressure sensor). In this study 5 alongshore arrays of pressure sensors (labled B, 
C, E, F, and 8m arrays, sensors indicated in blue) are used 
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Figure 3. Three-dimensional view of pressure sensor arrays used in this study. 
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Figure 4. Sonar Altimeter Locations. 
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Figure 5. Comparison of beach profiles obtained from CRAB and sonar 
altimeter data. The sonar measurements (asterisks) were interpolated linearly, except 
over the sand bar where a cubic spline was used. (a) small wave conditions. (b) storm 
conditions. 
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Figure 6. (a) Evolution of offshore Hs (measured in 8m depth) during the four 


month deployment. (b) Root-mean-square difference between beach profile based on 
CRAB and sonar altimeter data. Vertical dashed lines indicate times of beach profiles 
shown in figure 5. 
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Figure 7. Frequency and frequency-directional spectra for Case I, (August 7, 1997 


01:00 EST). The top panels display the (direction-integrated) frequency spectra at each 
array (blue line is model prediction, red line is observed spectrum at that array). The 
middle and lower panels contain the observed and predicted frequency-directional 
spectra, respectively. Note the exclusion of high frequency, large angle waves from the 
initializing model spectrum at the F array. 
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Figure 8. Case I, cross-shore shoaling evolution of significant wave height Hs, 


asymmetry, and skewness. Predictions of model runs initialized in 8m (dashed) and 5m 
(solid) depths are compared to observations (*) at each array. The depth profile and array 
locations are shown in bottom panel. | 
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Figure 9. Case I, cross-shore shoaling evolution of infragravity variance, mean 
direction, and directional spread. Predictions of model runs initialized in 8m (dashed) 
and 5m (solid) depths are compared to observations (*) at each array. The depth profile 
and array locations are shown in bottom panel. 
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Figure 10. Frequency and frequency-directional spectra for Case II, (August 10, 1997 
04:00 EST) (same format as figure 7). 
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Figure 11. Case II, cross-shore shoaling evolution of significant wave height Hs, 


asymmetry, and skewness (same format as figure 8). 
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Figure 12. Case II, cross-shore shoaling evolution of infragravity variance, mean 
direction, and directional spread (same format as figure 9). 
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Figure 13. Frequency and frequency directional spectra for Case III (September 12, 


1997 07:00 EST) (same format as figure 7). In order to enhance the directional properties 
at higher frequencies, the frequency-directional spectra (bottom two panels) were 
multiplied by the frequency. 
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Figure 14. Case III, cross-shore shoaling evolution of significant wave height Hs, 
asymmetry, and skewness (same format as figure 8). 
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Figure 15. Case III, cross-shore shoaling evolution of infragravity variance, mean 
direction, and directional spread (same format as figure 9). 
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Figure 16. Frequency and frequency-directional spectra for Case IV (October 19, 
1997 16:00 EST) (same format as figure 7). Note the exclusion of high frequency, large 
angle waves from the initializing model spectrum at the F array. In order to enhance the 

directional properties at higher frequencies, the frequency-directional spectra (bottom two 
panels) were multiplied by the frequency. 
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Figure 17. Case IV, cross-shore shoaling evolution of significant wave height Hs, 
asymmetry, and skewness (same format as figure 8). 
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Figure 18. Case IV: cross-shore shoaling evolution of infragravity variance, mean 
direction, and directional spread (same format as figure 9). 
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Figure 19. Left panels: predicted versus observed significant wave height (Hs) at the 
E, C, and B arrays. The model was initialized at the F array. Right panels: predicted 
versus observed significant wave height (Hs) both normalized by the initial Hs value. 

Solid diagonal line indicates a one-to-one correspondence, dashed line is a linear 
regression fit to all data points. 
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Figure 20. Predicted versus observed infragravity variance, both normalized by the 


initial infragravity variance. Left panels: results at the F, E, C, and B arrays for model 
runs initialized at the 8m array. Right panels: results at the E, C, and B arrays for model 
runs initialized at the F array. Solid diagonal line indicates a one-to-one correspondence, 
dashed line is a linear regression fit to all data points. 
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Figure 21. Predicted versus observed skewness (left panels) and asymmetry (right 


panels) at the E, C, and B arrays. The model was initialized at the F array. Solid diagonal 
line indicates a one-to-one correspondence, dashed line is a linear regression fit to all data 
points. 
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Figure 22. Predicted versus observed mean direction (left panels) and directional 


spread (right panels) at the E, C, and B arrays. The model was initialized at the F array. 
Solid diagonal line indicates a one-to-one correspondence, dashed line is a linear 
regression fit to all data points. 
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8m Model Run | Bias | = RMS-Error | ScatterIndex _| Correlation Coef._ 
| 0.092 | 0.986 

oo | oof | 043 | 0903 
01s =| aa | 7 9 
| tg | 0.976 

Mean Direction | —14deq_ | = =59deq | 04 | 0.915 
| zt | 0.792 











Table 1. Error statistics for model predictions initialized at the 8m array. Results 
are included from the F,E,C, and B arrays. 





F Model Run | Bias | — RMSError | ScatterIndex | Correlation Coef.| 


ee ey 

Asymmetry ss | —— 0.024 0.047 | 0.45 | 05" 
ISkewness sd] i005 S| S| 8 0.947 
-0.00034m42 | 0.00070m‘2_ | 021 &'||  o98 | 


a | 38deq | 0.28 | 0.972 


| 0972. 
Directional Spread | 3.4deq 3.8 dec 0.807 







Table 2. Error statistics for model predictions initialized at the F array. Results are 
included from the E,C, and B arrays. 
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